{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第四章　马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC）\n",
    "\n",
    "\n",
    "作者：[王何宇](http://person.zju.edu.cn/wangheyu)\n",
    "\n",
    "[浙江大学数学科学学院](http://www.math.zju.edu.cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们现在考虑如何为一个随时间变化的系统建模。首先还是需要一些基本假设。\n",
    "## 马尔可夫链\n",
    "假设一个随时间变化的系统，只有有限个状态，分别记作\n",
    "$$\n",
    "\\Omega = \\{x_1, x_2, \\cdots, x_n\\}.\n",
    "$$\n",
    "然后这里时间也是离散整点的，于是我们需要考虑的就是在 $t = 1, 2, \\cdots$ 时刻，系统的状态 $X_t \\in \\Omega$。如果，$X_t$ 的取值仅由 $X_{t - 1}$ 决定，和更早时刻的系统状态无关（**遗忘性**），则称满足此要求的状态序列 $\\{X_t\\}, t = 0, 1, \\cdots$ 为一个**马尔可夫链**。这里 $X_t$ 故意大写，是为了表示它是一个不确定的随机值，它可能是 $\\Omega$ 中任何一种状态，它实际上表示的是在 $\\Omega$ 中取值的一种随机分布，$\\Omega$ 中每一种状态，都以一个确定（仅在 $t$ 时刻）的概率取到。严格地说，马尔可夫链是一种随机分布序列，$t$ 时刻的分布 $X_t$ 仅由上一时刻 $X_{t - 1}$ 的分布决定。因此，如果我们确定给一个 $t - 1$ 时刻的系统状态，也即 $X_{t - 1}$ 的值，比如\n",
    "$$\n",
    "X_{t - 1} = x_i \\in \\Omega,\n",
    "$$\n",
    "那么，$X_t$ 的**分布**就只能有这个 $x_i$ 状态所**完全**确定。\n",
    "记\n",
    "$$\n",
    "P\\{x_i \\to x_j\\} := p_{ij}, \\quad i, j = 1, 2, \\cdots, n.\n",
    "$$\n",
    "则显然有\n",
    "$$\n",
    "\\sum_{j = 1}^n p_{ij} = 1, \\quad i = 1, 2, \\cdots, n.\n",
    "$$\n",
    "（这就是“完全”的意义。）\n",
    "\n",
    "形式上，这个矩阵\n",
    "$$\n",
    "P := \\left[p_{ij}\\right]_{n \\times n}, \\quad i, j = 1, 2, \\cdots, n.\n",
    "$$\n",
    "包含了一个马尔可夫链的全部信息，它的第 $i$ 行，代表了状态 $x_i$ 在下一个时刻可以转变的全部**可能**状态的分布，这里我们指出：\n",
    "1. 对角元 $p_{ii}$ 可以大于零，也即存在一定的概率，在下一个时间系统保持状态 $x_i$ 不变；\n",
    "2. 允许 $p_{ij} = 0$，也即存在这样的状态 $x_j$，它不可能由 $x_i$ 在一个时间步转变而来；\n",
    "3. $p_{ij} > 0$，$i, j = 1, 2, \\cdots, n$.\n",
    "\n",
    "注意这里我们实际上约定了 $p_{ij}$ 是一个和时间 $t$ 无关的常量。这样的马尔可夫链，称为齐次的（homogenous）。在接下去的讨论中，我们将模型局限在有限、齐次马尔可夫链下。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 令 $\\Omega = \\{x_1, x_2, x_3, x_4\\}$，转移矩阵为\n",
    "$$\n",
    "\\mathbf{P} = \\left[\n",
    "\\begin{array}{cccc}\n",
    "p_{11} & p_{12} & 0 & p_{14} \\\\\n",
    "p_{21} & 0 & p_{23} & 0 \\\\\n",
    "p_{31} & 0 & 0 & p_{34} \\\\\n",
    "0 & 0 & p_{43} & p_{44}\n",
    "\\end{array}\n",
    "\\right].\n",
    "$$\n",
    "或者等价地，可以用转移图表示为：\n",
    "<img src=\"graph.pdf\" width=\"25%\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们假定初始状态为\n",
    "$$\n",
    "X_0 = x_1,\n",
    "$$\n",
    "则初始状态是确定的，也即\n",
    "$$\n",
    "\\mathbf{p}_0 = (1, 0, 0, 0).\n",
    "$$\n",
    "而下一个状态 $X_1 = x_j$ 的概率分布是 $p_{1j}$，也即\n",
    "$$\n",
    "P\\{X_1 = x_j\\} = p_{1j}, \\quad j = 1, 2, 3, 4.\n",
    "$$\n",
    "\n",
    "一般地，若用\n",
    "$$\n",
    "\\mathbf{p}_t = (p_1, p_2, p_3, p_4)\n",
    "$$\n",
    "来表示 $t$ 时刻系统状态变量 $X_t$ 的概率质量分布（PMF），即\n",
    "$$\n",
    "\\mathbf{p}_j = P\\{X_t = x_j\\}, \\quad j = 1, 2, 3, 4,\n",
    "$$\n",
    "则称向量 $\\mathbf{p}$ 为概率向量（一般约定，它是个行向量）。由马尔可夫链定义，知\n",
    "$$\n",
    "\\mathbf{p}_{t + 1} = \\mathbf{p}_t \\mathbf{P}.\n",
    "$$\n",
    "注意这里行向量用法，进而\n",
    "$$\n",
    "\\mathbf{p}_{t + 1} = \\mathbf{p}_t\\mathbf{P} = \\mathbf{p}_{t - 1}\\mathbf{P}^2 = \\cdots = \\mathbf{p}_0 \\mathbf{P}^{t - 1}.\n",
    "$$\n",
    "到此我们可以给出结论：一个马尔可夫链的演化过程，完全由系统初始状态和转移矩阵确定。由上一章所学，我们可以在一个实际模拟中，通过随机抽取，来模拟一个马尔可夫链。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 假设股市是一个有限齐次马尔可夫链（韭菜田），并记三种状态分别为涨、平、跌对应 $x_1$, $x_2$ 和 $x_3$，其转移矩阵为\n",
    "$$\n",
    "\\mathbf{P} = \\left[\n",
    "\\begin{array}{ccc}\n",
    "0.3 & 0.2 & 0.5 \\\\\n",
    "0.4 & 0.2 & 0.4 \\\\\n",
    "0.4 & 0.3 & 0.3\n",
    "\\end{array}\n",
    "\\right],\n",
    "$$\n",
    "于是只要有一个初始状态，我们就可以对股市进行模拟。比如，记初始状态为涨，即\n",
    "$$\n",
    "\\mathbf{p}_0 = (1, 0, 0),\n",
    "$$\n",
    "则经过充分长时间的变化后，会发生什么？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "# 这句话是非标准的python，用于ipthon或jupyter这样的系统中，表示绘图即刻自动展开。\n",
    "%matplotlib inline\n",
    "\n",
    "# 这里把全部Warning过滤掉了. \n",
    "# 参见https://docs.python.org/2/library/warnings.html\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "#np.random.seed(250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们设置初值和转移矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = np.array([1, 0, 0])    \n",
    "P = np.array([[0.3, 0.2, 0.5], [0.4, 0.2, 0.4], [0.4, 0.3, 0.3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后每一次状态变换，都是一次向量乘以矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = p @ P\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "重复上述操作多次，最后我们发现，\n",
    "$$\n",
    "\\mathbf{p}_t \\to (0.36363636, 0.23966942, 0.39669421), \\quad t \\to \\infty.\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一般地，如果马尔可夫链链的概率向量 $\\mathbf{p}_t$ 当 $t \\to \\infty$ 时趋于一个常向量，我们称其为马尔可夫链的**不变分布**，或**平稳分布**，或**均衡分布**。自然界中有很多现象，是满足这个模型的，但是股市显然不是。\n",
    "\n",
    "从不变分布的定义，我们可以看到，事实上，不变分布就是迭代过程\n",
    "$$\n",
    "\\mathbf{p}_{t + 1} = \\mathbf{p}_t \\mathbf{P}\n",
    "$$\n",
    "的不动点 $\\pi$，满足\n",
    "$$\n",
    "\\pi \\mathbf{P} = \\pi,\n",
    "$$\n",
    "或者就是转移矩阵 $\\mathbf{P}$ 的特征值为 $1$ 对应的左特征向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理（Perron-Frobenius）** 若转移矩阵 $\\mathbf{P}$ 满足\n",
    "$$\n",
    "p_{ij} > \\delta > 0, \\quad i, j = 1, 2, \\cdots, n, \n",
    "$$\n",
    "则\n",
    "1. $\\mathbf{P}$ 存在特征值 $1$，且对应的左特征向量 $\\mathbf{w}$ 全部分量严格为正；\n",
    "2. 若 $\\mathbf{w}$ 归一化，则\n",
    "$$\n",
    "\\lim_{t \\to \\infty}\\mathbf{P}^t = \\mathbf{1w};\n",
    "$$\n",
    "3. 对于任何概率（行）向量 $\\alpha$，有\n",
    "$$\n",
    "\\lim_{t \\to \\infty} \\alpha \\mathbf{P}^t = \\mathbf{w},\n",
    "$$\n",
    "这里 $\\mathbf{w}$ 是已经归一的，$\\mathbf{1}$ 是分量均为 $1$ 的列向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理证明略，我们这里简单说明一下为何 $\\mathbf{P}$ 必然有特征值 $1$。事实上，令 $\\lambda$ 是 $\\mathbf{P}$ 的特征值，$\\mathbf{w} = (w_1, w_2, \\cdots, w_3)$ 是其左特征向量，即：\n",
    "$$\n",
    "\\mathbf{w} \\mathbf{P} = (\\sum_{i = 1}^nw_ip_{i1}, \\sum_{i = 1}^nw_ip_{i2}, \\cdots, \\sum_{i = 1}^nw_ip_{in}) = \\lambda \\mathbf{w} = (\\lambda w_1, \\lambda w_2, \\cdots, \\lambda w_n)\n",
    "$$\n",
    "将两边的各分量加在一起，并代入转移矩阵定义\n",
    "$$\n",
    "\\sum_{j = 1}^n p_{ij} = 1, \\quad i = 1, 2, \\cdots, n,\n",
    "$$\n",
    "得\n",
    "$$\n",
    "\\sum_{i = 1}^n w_i = \\lambda \\sum_{i = 1}^n w_i \\Rightarrow \\lambda = 1.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里注意，定理的要求是各元素均本质大于零，这是一个很强的条件。也就是任何两种状态，在一步之内，均有可能互变。实际上，这个条件可以放宽到：\n",
    "1. 系统的任何两个状态，都能在有限步互变，这一点又称为是**不可约**的；\n",
    "2. 系统的状态概率分布（注意不是状态）不存在周期变化；\n",
    "任何一个非周期不可约的马尔可夫链，都存在平稳分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之前的推导已经明示了平稳分布的求法：\n",
    "+ 理论上，它可以通过计算转移矩阵的 $1$ 特征值对应的左特征向量得到；\n",
    "+ 计算上，它可以通过不断抽样模拟马尔可夫链（迭代）得到。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC 抽样——Metropolis 算法\n",
    "\n",
    "我们继续抽样的话题。之前我们已经学习到：\n",
    "1. 如果存在累计分布函数的逆映射，则我们可以直接抽取指定分布；\n",
    "2. 如果只有概率密度函数，无法求得累计分布函数的逆映射，那么我们可以使用拒绝-接受法（AR）抽取指定分布；\n",
    "但是在一些实际工作中，我们只知道某个事件是随机的，并且知道它的一些分布规律，但是连概率密度函数 $f(x)$ 都无法给出。因为概率密度函数，必须满足\n",
    "$$\n",
    "\\int_{\\infty}^{\\infty} f(x) dx = 1,\n",
    "$$\n",
    "\n",
    "而很多时候，我们只知道某个随机事件，服从概率分布规律 $f(x)$，但是 $f(x)$ 无法写出显式表达式，或者有表达式，但是无法求积分，此时我们就没有办法对 $f(x)$ 做归一化，得到概率密度函数。或者说，在应用 AR 法时，本质困难出现在我们不知道 $f(x)$ 的上界，它甚至没有上界，就是一个一直增长的随机事件。\n",
    "\n",
    "在这种情况下，我们仍然可以利用马尔可夫链来获得指定随机事件的抽取（随机模拟）。基本想法是，我们构建一个马尔可夫链，使得它的平稳分布趋于指定的随机事件。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先，我们要假设，在指定的状态空间 $\\Omega = \\{x_1, x_2, \\cdots, x_n\\}$ 上，能够产生出我们的目标分布的概率密度函数，这里主要指状态空间要一致。然后关键在于，如果当前的状态是 $X_t$，那么下一个时刻的状态 $X_{t + 1}$ 如果确定，或者说，在已知当前状态为 $X_t$ 的前提下，$X_{t + 1}$ 满足什么概率分布。我们用概率论的记号，将这个分布记为 $g(\\cdot | X_t)$，称为**建议分布**。对应状态空间 $\\Omega$ 中的一个具体的状态 $x_i$, $i = 1, 2, \\cdots, n$，函数 $g$ 不必在整个 $\\Omega$ 上都本质大于零，实际上，$\\forall x \\in \\Omega$，我们称满足 $g(y | x) > 0$ 的状态构成的子集，为状态 $x$ 的**邻域**，记作 $\\mathcal{N}_x$。在很多算法中，我们要求 $g$ 满足对称性，也即\n",
    "$$\n",
    "g(y | x) = g(x | y), \\quad \\forall x, y \\in \\Omega.\n",
    "$$\n",
    "这一点在模型中往往会自然地做到（符合物理实际），并且它实际上可以弱化为\n",
    "$$\n",
    "g(y | x) > 0 \\Leftrightarrow g(x | y) > 0,  \\quad \\forall x, y \\in \\Omega,\n",
    "$$\n",
    "也即任两种状态的转移是可逆的。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去我们就可以逐一通过 $g$ 来模拟马尔可夫链的转移，但每一次从 $x$ 抽取到 $y$，我们都要做一个判定，到底需不需要转。也就是说，即便有 $X_t = x$，实际通过 $g(\\cdot | x)$ 抽取得到了 $y$，$X_{t + 1}$ 是否就取 $y$，还需要做一个独立抽取来判定，这里首先定义接受概率：\n",
    "$$\n",
    "h(x, y) = \\min \\left\\{1, \\frac{f(y)g(x | y)}{f(x)g(y | x)}\\right\\},\n",
    "$$\n",
    "特别地，当 $g$ 满足对称性条件时，上式简化为\n",
    "$$\n",
    "h(x, y) = \\min \\left\\{1, \\frac{f(y)}{f(x)}\\right\\}.\n",
    "$$\n",
    "我们首先计算这个概率，然后抽取一个 $(0, 1)$ 上的均匀分布的随机数 $r$，只有当 $r < h(x, y)$ 时，我们才令 $X_{t + 1} = y$，否则，仍然令 $X_{t + 1} = x$。也就是说，在这一个时间步，依概率 $h(x, y)$ 决定马尔可夫链是否从 $x$ 转移到 $y$。（这一个判定的意义在于，如果 $y$ 是一个概率极低的状态，那么它很难被转移到，但并不是完全没有希望。）\n",
    "这里，$f$ 只需要可计算就行，实际不需要满足概率密度函数归一化的要求（当然保正还是必须的）。而 $g$ 必须是一个严格的概率分布，但它是由我们设计的。这一个算法，被称为 Metropolis-Hastings 算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 我们先直接模拟两颗骰子的情形。显然两颗骰子的状态空间是 $\\Omega = \\{2, 3, \\cdots, 12\\}$。这个模型当然是可以得到概率质量分布的，但是我们这里直接通过 Metropolis-Hasting 算法直接产生符合这个模型的随机数。我们首先确定一下每一种状态可能的组合数，它实际上就对应这里的 $f(x)$，显然，\n",
    "$$\n",
    "f(2) = 1, f(3) = 2, f(4) = 3, f(5) = 4, f(6) = 5, f(7) = 6, f(8) = 5, f(9) = 4, f(10) = 3, f(11) = 2, f(12) = 1.\n",
    "$$\n",
    "注意这里 $f(x)$ 并未归一化。接下去考虑如何取邻域，也就是对 $x \\in \\Omega$，下一步 $y$ 可能的状态分布。这个取法可以很灵活，比如，取**极小邻域**，也就是 $\\forall x \\in \\Omega$，$\\mathcal{N}_x = \\{x - 1, x + 1\\}$，进一步加上对称性要求，我们就可以取\n",
    "$$\n",
    "g (y | x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\\frac{1}{2}, & y = \\max\\{x - 1, 2\\}, \\\\\n",
    "\\frac{1}{2}, & y = \\min\\{x + 1, 12\\},\\\\\n",
    "0, & \\mbox{others}.\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "这里注意考虑了 $\\Omega$ 两端的边界情况：\n",
    "$$\n",
    "g(3 | 2) = g(2 | 2) = g(11 / 12) = g(12 | 12) = \\frac{1}{2}. \n",
    "$$\n",
    "保证了对称性。这个自己设计的概率密度函数（严格说应该是概率质量分布）也可以等价地用下面的矩阵表示：\n",
    "$$\n",
    "G = \\left(\n",
    "\\begin{array}{ccccccc}\n",
    "\\frac{1}{2} & \\frac{1}{2} & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "\\frac{1}{2} & 0 & \\frac{1}{2} & \\cdots & 0 & 0 & 0 \\\\\n",
    "0 & \\frac{1}{2} & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots & \\vdots \\\\\n",
    "0 & 0 & 0 & \\cdots & 0 & \\frac{1}{2} & 0 \\\\\n",
    "0 & 0 & 0 & \\cdots & \\frac{1}{2} & 0 & \\frac{1}{2} \\\\\n",
    "0 & 0 & 0 & \\cdots & 0 & \\frac{1}{2} & \\frac{1}{2}\n",
    "\\end{array}\n",
    "\\right).\n",
    "$$\n",
    "这里分布的对称性和矩阵的对称性是一致的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20000\n",
    "f = np.array([0, 0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1])  \n",
    "d = np.zeros(N)\n",
    "x = 5\n",
    "for i in range(N):\n",
    "    U = np.random.rand()\n",
    "    if x is 2:\n",
    "        if U < 0.5:\n",
    "            y = 3\n",
    "        else:\n",
    "            y = 2\n",
    "    elif x is 12:\n",
    "        if U < 0.5:\n",
    "            y = 11\n",
    "        else:\n",
    "            y = 12\n",
    "    else:\n",
    "        if U < 0.5:\n",
    "            y = x - 1\n",
    "        else:\n",
    "            y = x + 1\n",
    "    h = np.min([1, f[y] / f[x]])\n",
    "    U = np.random.rand()\n",
    "    if U < h:\n",
    "        x = y;\n",
    "    d[i] = x        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(d, bins = 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "符合双骰子线性增长的规律。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "和极小邻域法类似，还有一种几乎万能的邻域取法是**极大邻域**法。这里我们直接取\n",
    "$$\n",
    "G = \\left(\n",
    "\\begin{array}{cccc}\n",
    "\\frac{1}{11} & \\frac{1}{11} & \\cdots & \\frac{1}{11} \\\\\n",
    "\\frac{1}{11} & \\frac{1}{11} & \\cdots & \\frac{1}{11} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\frac{1}{11} & \\frac{1}{11} & \\cdots & \\frac{1}{11} \\\\\n",
    "\\end{array}\n",
    "\\right),\n",
    "$$\n",
    "这里实际上取 $\\mathcal{N}_x = \\Omega$, 而 $g(\\cdot | x)$ 取了一个常量，这样的函数当然是对成的。这个模拟实际上写起来更简单："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 5\n",
    "for i in range(N):\n",
    "    y = np.random.randint(2, 13)\n",
    "    h = np.min([1, f[y] / f[x]])\n",
    "    U = np.random.rand()\n",
    "    if U < h:\n",
    "        x = y;\n",
    "    d[i] = x        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(d, bins = 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个算法，也可以使用在连续分布抽样上。比如，已知密度函数为\n",
    "$$\n",
    "f(x) = 0.5 x^2 e^{-x}\n",
    "$$\n",
    "的伽马分布，这里我们不用考察这个密度函数是否归一，它的状态空间是 $\\Omega = [0, +\\infty)$，所以我们在产生建议分布 $g$ 时也要注意，$g(y | x)$ 必须确保 $y > 0$，同时我们也希望继续保持 $g$ 的对称性。为此，我们取 $\\mathcal{N}_x = (x - 1, x + 1)$, 且 $g(\\cdot | x)$ 为其上的均匀分布。如果这里 $y < 0$，则令 $y = x$，这个做法类似之前离散的情况。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 0.5 * x * x * np.exp(-x) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 40000\n",
    "d = np.zeros(N)\n",
    "x = 2\n",
    "for i in range(N):\n",
    "    y = np.random.rand() * 2 - 1 + x\n",
    "    if y < 0:\n",
    "        y = x\n",
    "    h = np.min([1, f(y) / f(x)])\n",
    "    U = np.random.rand()\n",
    "    if U < h:\n",
    "        x = y\n",
    "    d[i] = x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num = plt.hist(d[20000:], bins = 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对比一下目标曲线，下图中的蓝线是根据统计结果做的归一化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = num[0]/(0.5 * N)\n",
    "x = num[1][0:-1] + np.diff(num[1]) * 0.5\n",
    "plt.plot(x,y)\n",
    "plt.hold\n",
    "plt.plot(x, f(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metropolis-Hasting 算法的有效性分析\n",
    "上面的例子展示了 Metropolis-Hasting 算法的有效性。它本质上可以看作是拒绝-接受法的一种扩展。但是对它为何会恰好正确，可能大家还是会有一些疑惑。我们继续用一种直观的方式解释这个现象。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们只考虑有限状态空间，在很多物理实际中，有限已经足够。比如我们考虑一个指定有限区域空间的气体密度变化，我们将这个空间打成细密均匀的网格，每一个网格中的分子数目，本质上就是气体的密度，同时可以看作一种概率状态分布。这是一个有限状态空间。我们用 $m_t(x)$ 表示 $t$ 时刻，空间 $x$ 位置（可以认为是一个空间坐标，或者一个网格编号）的分子数。那么到 $t + 1$ 时刻，$x$ 格子和一个和它相邻的 $y$ 格子之间，会有分子交换，导致 $x$ 格子中气体分子数的变化（在这一过程中，我们只考虑由 $x$ 和 $y$ 两个格子之间交换分子引起的数量变化），这一变化记为：\n",
    "$$\n",
    "\\delta_{xy} = m_t(x) P\\{x \\to y\\} - m_t(y) P\\{y \\to x\\} = m_t(y)P\\{x \\to y\\}\\left(\\frac{m_t(x)}{m_t(y)} - \\frac{P\\{y \\to x\\}}{P\\{x \\to y\\}}\\right).\n",
    "$$\n",
    "这里只用到了质量守恒定律。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据热力学第二定律，在无外界干扰的情况下，经过充分长的时间，气体状态会达到平衡状态，也就是从概率意义上，$\\delta_{xy} = 0$，也即 $\\forall x \\in \\Omega$，有\n",
    "$$\n",
    "\\delta_{xy} = 0 \\Rightarrow \\lim_{t \\to \\infty} \\frac{m_t(x)}{m_t(y)} = \\frac{P\\{y \\to x\\}}{P\\{x \\to y\\}},\n",
    "$$\n",
    "记\n",
    "$$\n",
    "m_e(x) = \\lim_{t \\to \\infty} m_t(x),\n",
    "$$\n",
    "则 $m_e$ 即气体在平衡态下的分布。由 Metropolis-Hasting 算法，\n",
    "$$\n",
    "P\\{x \\to y\\} = g(y | x) h(x, y),\n",
    "$$\n",
    "再由 $g$ 的对称性，有\n",
    "$$\n",
    "\\frac{m_e(x)}{m_e(y)} = \\lim_{t \\to \\infty} \\frac{m_t(x)}{m_t(y)} = \\frac{P\\{y \\to x\\}}{P\\{x \\to y\\}} = \\frac{h(y, x)}{h(x, y)},\n",
    "$$\n",
    "这里最后一项约去了 $g(y | x) = g(x | y)$。而对 $f(x) > f(y)$，$h(y, x) = 1$，此时\n",
    "$$\n",
    "\\frac{m_e(x)}{m_e(y)} = \\frac{1}{h(x, y)} = \\frac{f(x)}{f(y)},\n",
    "$$\n",
    "而对 $f(x) \\leq f(y)$，$h(x, y) = 1$，有\n",
    "$$\n",
    "\\frac{m_e(x)}{m_e(y)} = h(y, x) = \\frac{f(x)}{f(y)},\n",
    "$$\n",
    "由 $x$ 和 $y$ 的任意性，必有 $m_e(x)$ 和 $f(x)$ 成正比例。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "事实上，对平衡态的另一种微观等价描述是：\n",
    "$$\n",
    "f(x)P\\{x \\to y\\} = f(y)P\\{y \\to x\\}, \\quad \\forall x, y \\in \\Omega.\n",
    "$$\n",
    "而 Metropolis-Hasting 算法满足该微观等价性：\n",
    "$$\n",
    "f(x)h(x, y) = f(y)h(y, x).\n",
    "$$\n",
    "这个其实是算法的设计初衷。最后，Metropolis-Hasting 算法产生的马尔可夫链满足：\n",
    "1. 不可约：由 $g$ 的定义知，状态转移可以遍历状态空间 $\\Omega$；\n",
    "2. 非周期：$f$ 至少存在一个最大值 $y^*$，满足其停留的概率是正概率，也即存在可能状态永远停留在 $y^*$。\n",
    "所以该算法产生的马尔可夫链，必存在平稳分布。\n",
    "\n",
    "以上的分析过程并不完全数学上严格，但给出了明确的直观思路。有助于我们理解和掌握此算法，以后辅以概率论技术，也不难严格化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理**  若 $X_t$ 是一个定义在有限状态空间 $\\Omega$ 上的不可约非周期马尔可夫链，其平稳分布为 $\\pi$，设 $\\xi : \\Omega \\to \\mathbb{R}$ 为任意映射，则\n",
    "$$\n",
    "E_\\pi[\\xi] = \\lim_{N \\to \\infty} \\frac{1}{N} \\sum_{i = 1}^N \\xi(x_i)\n",
    "$$\n",
    "由这个定理，我们可以利用 MCMC 方法计算各种随机模型的期望。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 统计力学\n",
    "\n",
    "大家可能已经注意到，计算机模拟领域的开拓者们，大部分是物理学家。实际上，计算机模拟这个研究方向的诞生，有三个基本条件：计算机的出现；概率论的成熟；以及物理学（核平）研究的需要。\n",
    "\n",
    "十九世纪以后，人类对物质本质的了解，已经到了原子内部。但即便如此，我们仍然无法直接从分子运动去预测宏观气体的现象，如密度、压力和温度的变化等等。但之前在 MCMC 的理论来源时，我们已经知道，所谓气体的宏观密度，从微观角度看，就是气体分子的空间分布。由于分子数目巨大，我们只要知道其概率分布，就可以对宏观气体密度作出足够好的预测和估算。在类似的思想下，Boltzmann 等人开创了统计力学，建立了一系列描述微观粒子概率分布的方程，比如著名的 Boltzmann 方程：\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t} f(x, v, t) + v \\cdot \\nabla_x f(x, v, t) = Q(f, f) (x, v, t).\n",
    "$$\n",
    "理论上说，对这些方程的统计模拟，可以完全代替宏观模型下的流体方程，如 Euler、Navier-Stikes 等流体力学方程，对流体的宏观性质作出预测和估计。同样的思想，还可以用在电、磁甚至量子力学领域。那么如何正确理解统计物理学模型，如何结合概率论直观给出正确且高效的抽样，并在现代计算机上实现，就称为一项现代科研领域不可缺少的手段。\n",
    "\n",
    "我们这里考虑一个最简单的统计物理模型——**Ising** 模型。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1920 年，德国人 Ernst Ising 在他的博士论文中第一次用统计力学模型解释了铁磁体磁性相变的物理现象。也就是铁磁性物质，存在一个被称为 “居里温度” 的相变点，当温度超过居里温度后，磁性迅速消失，降温至居里温度之下，磁性又会恢复。Ising（实际上已经被我们大大简化）用一个空间网格来模拟铁磁性物质，每个网格代表一个微观铁原子的磁极状态，这个状态这里我们只取两种：北极向上记作 $1$，南极向上记作 $-1$。若大部分铁原子具有相同的磁极状态，那么物质整体就呈现宏观磁性；若原子的磁极是均匀混乱分布的，那么微观磁极互相抵消，宏观物质整体就没有磁性。这里我们假设原子总数是 $N$，那么不同原子的不同磁性选择，会形成 \n",
    "$$\n",
    "K = 2^N\n",
    "$$ 种可能的状态，这是一个巨大的数字，但仍然是有限的。我们这里给出假设：**每一种微观组合状态，都是等概出现的。** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了描述一个系统的混乱程度，Boltzmann 引入了熵（entrpy）的概念：\n",
    "$$\n",
    "S = k_B \\ln K,\n",
    "$$\n",
    "这里 $k_B$ 称为 Boltmann 常数。由热力学第二定律我们知道，熵总是增加的，因此一个平衡系统，熵是极大的（不能再乱了）。我们现在考虑有两个独立的系统，各自的状态数分别是 $K_1$ 和 $K_2$，现在观念上把两个系统合并成一个系统考虑，但实际物理上不考虑它们之间由分子交换引起的新状态，则总系统的状态数为：\n",
    "$$\n",
    "K = K_1 \\cdot K_2,\n",
    "$$\n",
    "而对应地，\n",
    "$$\n",
    "S = k_B \\ln(K_1 \\cdot K_2) = k_B \\ln K_1 + k_B \\ln K_2 = S_1 + S_2,\n",
    "$$\n",
    "即熵有可加性。再由能量守恒定律：\n",
    "$$\n",
    "E = E_1 + E_2 \\Rightarrow E_2 = E - E_1.\n",
    "$$\n",
    "现在假设总能量 $E$ 不随时间变化（但 $E_1$ 和 $E_2$ 会有彼增此减的变化），并将 $S$ 看作 $E_1$ 的函数（也就是看 $E_1$ 如何变化，使得 $S$ 达到最大），则有\n",
    "$$\n",
    "\\left.\n",
    "\\begin{array}{r}\n",
    "\\displaystyle 0 = \\frac{d S}{d E_1} = \\frac{d S_1}{d E_1} + \\frac{d S_2}{d E_2} \\frac{d E_2}{d E_1} \\\\\\\\\n",
    "\\displaystyle \\frac{d E_2}{d E_1} = \\frac{d E}{d E_1} - \\frac{d E_1}{d E_1} = -1\n",
    "\\end{array}\n",
    "\\right\\} \\Rightarrow \\frac{d S_1}{d E_1} - \\frac{d S_2}{d E_2} = 0 \\Rightarrow \\frac{d S_1}{d E_1} = \\frac{d S_2}{d E_2} \\left( = \\frac{d S}{d E}\\right)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过上面的推导，我们的到一个结论：均衡态的 $\\frac{d S}{d E}$ 在系统处处是常数，我们将它记作温度 $T$：\n",
    "$$\n",
    "\\frac{1}{T} := \\frac{d S}{d E}.\n",
    "$$\n",
    "热力学第二定律的另一种表达就是能量总是从高温系统流向低温系统，直至恒温。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去要考虑，对于一个处于平衡态的系统，其中微粒的实际状态服从什么分布。这里首先需要一个基本假设：**只要时间足够长，平衡系统中会经历全部可能的微观状态。**现在我们考虑由 $N$ 个互相不交换粒子的独立系统构成的**正则系综（canonical ensemble）**。它们的总能量是各个系统之和\n",
    "$$\n",
    "E = \\sum_{j = 1}^N E_j,\n",
    "$$\n",
    "总微观状态数是各系统微观状态数之积\n",
    "$$\n",
    "K = \\prod_{j = 1}^N K_j.\n",
    "$$\n",
    "我们令各系统都是保持恒温 $T$ 的，并且系统可能的微观状态共有 $m$ 个，分别是 $\\{\\psi_1, \\psi_2, \\cdots, \\psi_m\\}$，则我们希望能够得到系统处于 $\\psi_i$ 上的概率，也就是概率质量分布。令 $n_i$ 表示处于 $\\psi_i$ 状态的系统数，$E_i$ 表示 $\\psi_i$ 系统所具有的能量，则有\n",
    "$$\n",
    "\\sum_{i = 1}^m n_i = N, \\quad \\sum_{i = 1}^m n_i E_i = E.\n",
    "$$\n",
    "而系综的总微观状态数为\n",
    "$$\n",
    "K = \\frac{N!}{n_1! n_2! \\cdots n_m!}.\n",
    "$$\n",
    "这里分子是系统总数的全排列，分母表示在这个排列下，处于同一状态下的系统应该视作一种情况，因为他们之间交换粒子不产生新的状态。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在由最大熵原理，求 $\\{n_i\\}$ 的组合形式，使 $S$ 最大，即\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{rcl}\n",
    "\\displaystyle \\max_{\\{n_i\\}}S & = & k_B \\ln K\\\\\n",
    "\\mathrm{s. t.} && \\sum_{i = 1}^m n_i = N, \\\\\n",
    "&& \\sum_{i = 1}^m n_iE_i = E.\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "由 Lagrange 乘子法，\n",
    "$$\n",
    "\\mathcal{L}\\left(\\{n_i\\}, \\alpha, \\beta\\right) = k_B \\ln K - \\alpha \\sum_{i = 1}^m n_i - \\beta \\sum_{i = 1}^m n_i E_i,\n",
    "$$\n",
    "其中 $\\alpha$ 和 $\\beta$ 是对应的 Lagrange 乘子。于是\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial n_i} & = & k_B\\frac{\\partial}{\\partial n_i} \\ln K - \\alpha \\frac{\\partial}{\\partial n_i} \\sum_{j = 1}^m n_j - \\beta \\frac{\\partial}{\\partial n_i}\\sum_{j = 1}^m n_j E_i \\\\\n",
    "& = & k_B \\frac{\\partial}{\\partial n_i} \\ln K - \\alpha - \\beta E_i = 0.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "注意这里第一行的第二项当 $i = j$ 时是 $1$，其余都是 $0$；第三项当 $i = j$ 时是 $E_i$，其余是 $0$。而对第二行第一项，由 Stirling 公式，\n",
    "$$\n",
    "\\ln n! \\approx n \\ln n - n,\n",
    "$$\n",
    "故\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\ln \\frac{N!}{n_1!\\cdots n_m!} & = & \\ln N! - \\left(\\sum_{i = 1}^m \\ln n_i\\right) \\\\\n",
    "& = & N \\ln N - N - \\sum_{i} n_i \\ln n_i + \\sum_{i} n_i\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "注意到 \n",
    "$$\n",
    "N = \\sum_i n_i,\n",
    "$$\n",
    "于是\n",
    "$$\n",
    "\\ln \\frac{N!}{n_1!\\cdots n_m!} = N \\ln N - \\sum_i n_i \\ln n_i.\n",
    "$$\n",
    "进而\n",
    "$$\n",
    "\\frac{\\partial \\ln K}{\\partial n_i} = \\frac{\\partial}{\\partial n_i}\\left(N\\ln N - \\sum_j n_j \\ln n_j\\right),\n",
    "$$\n",
    "这里第一项是常数，求导后为零；第二项根据求导乘法展开，只有 $i = j$ 才非零，故\n",
    "$$\n",
    "\\frac{\\partial \\ln K}{\\partial n_i} = -\\frac{\\partial}{\\partial n_i} n_i \\ln n_i = -\\ln n_i - 1,\n",
    "$$\n",
    "综合有，\n",
    "$$\n",
    "-k_B \\ln n_i - k_B - \\alpha - \\beta E_i = 0,\n",
    "$$\n",
    "所以\n",
    "$$\n",
    "n_i = e^{-1 - \\frac{\\alpha}{k_B} - \\beta \\frac{E_i}{k_B}},\n",
    "$$\n",
    "和\n",
    "$$\n",
    "e^{-\\frac{\\beta}{k_B}E_i}\n",
    "$$\n",
    "成正比。因此，系统在各状态 $\\Omega = \\{\\psi_1, \\psi_2, \\cdots, \\psi_N\\}$ 上的概率密度函数（这里其实还是概率质量）为\n",
    "$$\n",
    "p_i = \\frac{e^{-\\frac{\\beta}{k_B}}E_i}{Z},\n",
    "$$\n",
    "其中 $Z$ 是归一化因子，满足\n",
    "$$\n",
    "Z = \\sum_{i = 1}^m e^{-\\frac{\\beta}{k_B}}E_i,\n",
    "$$\n",
    "使得\n",
    "$$\n",
    "\\sum_{i \\in \\Omega}p_i = 1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去我们估计 $\\beta$，为此我们首先修改一下模型中熵的定义：\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{rcl}\n",
    "\\displaystyle \\max_{n_i} S & = & -k_B \\sum_{i = 1}^m p_i \\log p_i \\\\\n",
    "\\mathrm{s. t.} && \\sum_{i = 1}^m p_i = 1\\\\\n",
    "&& \\sum_{i = 1}^m p_i E_i = \\bar{E}\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "这里采用的熵定义是信息论大佬 Shannon 引入的信息熵：\n",
    "$$\n",
    "S = - \\sum_{i = 1}^m p_i \\ln p_i,\n",
    "$$\n",
    "$p_i$ 是某概率分布。它和 Boltzmann 引入的熵是等价的（差一个常数因子）。我们仍然使用 Lagrange 乘子法来分析：\n",
    "$$\n",
    "\\mathcal{L}\\left(\\{p_i\\}, \\alpha, \\beta\\right) = -k_b \\sum_{i = 1}^m p_i \\ln p_i - \\alpha \\sum_{i = 1}^m p_i - \\beta \\sum_{i = 1}^m p_i E_i,\n",
    "$$\n",
    "则\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial p_i} = \\sum_{i = 1}^m (-k_B \\ln p_i - k_B - \\alpha - \\beta E_i) = 0,\n",
    "$$\n",
    "对任意的 $\\{p_i\\}$ 均成立，因此\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "&& -k_B \\ln p_i - k_B - \\alpha - \\beta E_i  =  0, \\\\\n",
    "&\\Rightarrow& -k_B \\ln p_i  =  k_B + \\alpha + \\beta E_i, \\\\\n",
    "&\\Rightarrow& -k_B p_i \\ln p_i  =  k_B p_i + \\alpha p_i + \\beta E_i p_i, \\\\\n",
    "&\\Rightarrow& -k_B \\sum_{i = 1}^m p_i \\ln p_i  =  \\sum_{i = 1}^m (k_B p_i + \\alpha p_i + \\beta E_i p_i), \\\\\n",
    "&& =  k_B + \\alpha + \\beta \\bar{E}, \\\\\n",
    "&\\Rightarrow& S = k_B + \\alpha + \\beta \\bar{E}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "由温度定义，在平衡态，由\n",
    "$$\n",
    "\\beta = \\frac{\\partial S}{\\partial \\bar{E}} = \\frac{1}{T}.\n",
    "$$\n",
    "代入 $p_i$ 的表达式，有\n",
    "$$\n",
    "p_i = \\frac{1}{Z}e^{-\\frac{E_i}{k_BT}}, \\quad Z = \\sum_{i = 1}^m e^{-\\frac{E_i}{k_BT}}.\n",
    "$$\n",
    "这个分布，被称为 Boltzmann 分布。它为我们微观模拟提供了抽样依据。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们终于可以回到 Ising 模型。铁磁体的总磁矩量是一个宏观量，但它事实上是由单个原子自旋产生的微观磁场的叠加。现在假设目标铁磁体处于热平衡状态，而各原子的微观量服从 Boltzmann 分布，而宏观量本质上就是这些微观量的一种加权叠加，或称期望：\n",
    "$$\n",
    "<A> = \\sum_{x \\in \\Omega} A_x\\frac{e^{-\\frac{E_x}{k_BT}}}{Z},\n",
    "$$\n",
    "其中 $A_x$ 代表在 $x$ 位置的微观量，$<A>$ 称为 Boltzmann 期望，对应我们观念中的宏观量。在有了 Metropolis 算法之后，我们可以不用纠结如何求归一化系数 $Z$，于是一切都可以实际抽样和模拟。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实际算例\n",
    "\n",
    "我们来考虑一个简化版的 Ising 模型。用一个 $N \\times N$ 的等距二维网格来模拟一块铁磁体。每个网格交叉点代表一个铁原子，它具有 $\\pm 1$ 两种自选状态之一。即\n",
    "$$\n",
    "s(i, j) = \\pm 1, \\quad \\Omega = \\{S | S = s(i, j),\\quad i, j = 1, 2, \\cdots, N\\},\n",
    "$$\n",
    "而总磁矩为\n",
    "$$\n",
    "M(S) = \\frac{1}{N^2}\\sum_{j = 1}^N \\sum_{i = 1}^N s(i, j).\n",
    "$$\n",
    "显然，若 $s(i, j)$ 是均匀分布的，那么 $M(S) \\approx 0$.\n",
    "\n",
    "对于一个实际模型，其边界，也就是铁磁体如何和外界交换信息是重要的。在本问题中，我们试图回避这个问题，或者，我们的理想情况就是一块无限边界的铁磁体，这样我们可以对有限区域的铁磁体采用**周期边界条件**：\n",
    "$$\n",
    "\\begin{array}{rclrcl}\n",
    "s(0, j) &=& s(N, j), & s(N + 1, j) &=& s(1, j),\\\\\n",
    "s(i, 0) &=& s(i, N), & s(i, N + 1) &=& s(i, 1),\\\\\n",
    "s(0, 0) &=& s(N, N), & s(0, N + 1) &=& s(N, 1),\\\\\n",
    "s(N + 1, 0) &=& s(1, N), & s(N + 1, N + 1) &=& s(1, 1).\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "这样相当于模拟一块无限但是结构循环出现的物质。我们也可以假设不存在外界，或者物质与外界彻底隔绝，表现为我们关心的量，比如这里的 $M$，延边界法向的变化率为零，也即\n",
    "$$\n",
    "\\frac{\\partial M}{\\partial \\vec{n}} = 0.\n",
    "$$\n",
    "这种边界也称为**自由边界条件**。在实际抽样时，直接当外部格点不存在就行了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在 $\\forall (i, j)$，定义邻域为\n",
    "$$\n",
    "\\Lambda (i, j) = \\{(i - 1, j), (i, j - 1), (i, j + 1), (i + 1, j)\\},\n",
    "$$\n",
    "也就是和 $(i, j)$ 直接相邻的格点。而能量则由经验公式\n",
    "$$\n",
    "E(S) = -\\frac{1}{2} J \\sum_{j = 1}^N \\sum_{i = 1}^N \\sum_{\\mu \\in \\Lambda (i, j)} s(\\mu) s(i, j),\n",
    "$$\n",
    "决定，这里 $J > 0$ 是常数。而由之前分析，当系统处于热平衡态时，状态为 $s \\in \\Omega$ 的概率服从 Boltzmann 分布，即\n",
    "$$\n",
    "P(S) = \\frac{1}{Z}e^{-\\frac{E(S)}{k_BT}}, \\quad S \\in \\Omega.\n",
    "$$\n",
    "因此我们不难用随机模拟来估计总磁矩量\n",
    "$$\n",
    "<M> = \\frac{1}{Z} \\sum_{S \\in \\Omega} M(S) e^{-\\beta E(S)}, \\quad \\beta = \\frac{1}{k_BT}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们用 Metropolis-Hasting 算法来模拟抽样，这里可以采用最小邻域（变化）思想，即每一次的建议改变，都是在现有的一个状态 $S$ 中随机抽取一个格点，然后对其微磁矩进行反转，得到一个新的状态 $S'$，然后根据\n",
    "$$\n",
    "h(S, S') = \\min \\left\\{1, \\frac{P(S')}{P(S)}\\right\\} = \\min\\left\\{1, \\frac{e^{-\\beta E(S')}/Z}{e^{-\\beta E(SS)}/Z}\\right\\} = \\min \\left\\{1, e^{-\\beta[E(S')-E(S)]}\\right\\}\n",
    "$$\n",
    "来判定是否接受建议改变。下面是模拟的代码。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 3\n",
    "S = np.random.randint(0, 2, (N, N)) * 2 - 1\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Energy(S):\n",
    "    E = S[0][0] * S[0][1] + S[0][0] * S[1][0]\n",
    "    E = E + S[N - 1][0] * S[N - 2][0] + S[N - 1][0] * S[N - 1][1]\n",
    "    E = E + S[0][N - 1] * S[0][N - 2] + S[0][N - 1] * S[1][N - 1]\n",
    "    E = E + S[N - 1][N - 1] * S[N - 1][N - 2] + S[N - 1][N - 1] * S[N - 2][N - 1]\n",
    "    for i in range(1, N - 1):\n",
    "        E = E + S[i][0] * S[i - 1][0] + S[i][0] * S[i + 1][0] + S[i][0] * S[i][1]\n",
    "        E = E + S[i][N - 1] * S[i - 1][N - 1] + S[i][N - 1] * S[i + 1][N - 1] + S[i][N - 1] * S[i][N - 2]\n",
    "        for j in range(1, N - 1):\n",
    "            E = E + S[0][j] * S[0][j - 1] + S[0][j] * S[0][j + 1] + S[0][j] * S[1][j]\n",
    "            E = E + S[N - 1][j] * S[N - 1][j - 1] + S[N - 1][j] * S[N - 1][j + 1] + S[N - 1][j] * S[N - 2][j]\n",
    "            E = E + S[i][j] * S[i - 1][j] + S[i][j] * S[i + 1][j] + S[i][j] * S[i][j + 1] + S[i][j] * S[i][j - 1]\n",
    "    E = -E * 0.5\n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Sample(warmup, trials, S, T):\n",
    "    M = np.zeros(trials)\n",
    "    beta = 1 / T\n",
    "    Es = Energy(S)\n",
    "    for t in range(warmup + trials):\n",
    "        i = np.random.randint(N)\n",
    "        j = np.random.randint(N)\n",
    "        Y = S.copy()\n",
    "        Y[i][j] = -Y[i][j]\n",
    "        Ey = Energy(Y)\n",
    "        H = np.min([1, np.exp(-beta * (Ey - Es))])\n",
    "        if np.random.rand() < H:\n",
    "            S = Y.copy()\n",
    "            Es = Ey\n",
    "        if t >= warmup:\n",
    "            M[t - warmup] = np.sum(S)\n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = Sample(100000, 100000, S, 1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(M, bins = 18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里有一个坑，Python 默认的数组复制是传址的，比如："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones(4)\n",
    "print(a)\n",
    "b = a\n",
    "b[1] = 9\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意到这里 a 和 b 一起被修改。正确的复制操作是对 list："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = [1, 2, 3, 4]\n",
    "b = a[:]\n",
    "b[2] = 9\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对 array："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3, 4])\n",
    "b = a.copy()\n",
    "b[2] = 9\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ising(T):    # T is the tempure.  +-+-+   Ising area.\n",
    "    nTrials = 100000    # for stat.          |3|4|   N = 2.\n",
    "    startup = 100000    # for steady.        +-+-+   free boundary.\n",
    "    d = np.zeros(nTrials)   # results.       |1|2| \n",
    "    beta = 1/T             #                +-+-+ \n",
    "    M = 0\n",
    "    s = 2 * np.random.randint(0, 2, 5) - 1   # randomly set 1 or -1.\n",
    "    Es = s[1] * s[2] + s[3] * s[4]        \n",
    "    Es = Es + s[1] * s[3] + s[2] * s[4]   \n",
    "    Es = -Es                              \n",
    "    for t in range(startup + nTrials):        \n",
    "        k = np.random.randint(1, 5)  # randomly pick 1 ~ 4.\n",
    "        #print(\"k, S = \", k, s, Es)\n",
    "        y = s.copy()         # y is the suggestion dist in s' neighbour. \n",
    "        y[k] = -s[k]  # randomly flip the status on one point.\n",
    "        #print(\"s = \",  s)\n",
    "        Ey = y[1] * y[2] + y[3] * y[4]\n",
    "        Ey = Ey + y[1] * y[3] + y[2] * y[4]\n",
    "        Ey = -Ey\n",
    "        #print(\"k, Y = \", k, y, Ey)\n",
    "        h = np.min([1, np.exp(-beta * (Ey - Es))])\n",
    "        if (np.random.rand() < h):\n",
    "            s = y[:]\n",
    "            Es = Ey\n",
    "            #print(s, Es)\n",
    "        if (t >= startup):\n",
    "            Ms = s[1] + s[2] + s[3] + s[4]  \n",
    "            M = M + Ms                      \n",
    "            d[t - startup] = Ms\n",
    "        \n",
    "    \n",
    "#    x = -4:1:4;\n",
    "    plt.hist(d);\n",
    "#    M = (M / nTrials) / 4;\n",
    "#    D = sum(d.^2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ising(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
